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Abstract 



These are lecture notes, in three main parts: phase space diagrams, Fermi-Eyges theory, and the 
application of the theory to transverse beam spreading in homogeneous matter. 

Phase space diagrams serve to visualize correlations between position and direction in a particle 
beam. They do not in themselves yield quantitative formulas, but they facilitate geometric reasoning 
which makes such formulas more plausible. We introduce phase space via a model beam line where 
multiple Coulomb scattering (MCS) and drift (motion along the beam line) happen separately 
(usually they are combined). We discuss the effect of a pure scatterer, and of a pure drift, on the 
phase space distribution. We introduce the beam ellipse and the area enclosed by it, the cmittance. 
We discuss how emittance changes in a drift and in a thin scatterer. Finally, we discuss phase space 
diagrams of a more realistic beam line. 

In a brief detour, we review various mathematical topics, justify the Gaussian approximation by 
comparing it to the full theory of MCS, and briefly discuss the concept 'scattering power' which is 
needed in Fermi-Eyges theory. 

When the phase space distribution happens to be Gaussian it can be computed by means of 
Fermi-Eyges theory. The most general case is a Gaussian incident beam traversing one or more 
homogeneous slabs of varying materials. All the resulting phase space distributions will be Gaussian 
until the beam stops or encounters a beam limiting device or a transverse heterogeneity. We discuss 
Eyges' original theory for an ideal beam entering a single slab, and its later generalization to a finite 
(but still Gaussian) beam entering a stack of slabs. We revisit the beam ellipse, showing it to be 
an iso-density contour in phase space. We use the theory to confirm the emittance change in a drift 
and in a thin scatterer. We discuss equivalent sources and finally, we compute the fraction of the 
beam contained by the beam ellipse defined in a more general way. 

Preston and Koehler (PK), without using Fermi-Eyges theory, studied transverse spreading of 
particle beams stopping in matter. We recast their derivations in Fermi-Eyges language. First, we 
show that their basic equation for the rms beam size at any measuring plane, which they obtained by 
direct physical reasoning, is equivalent to a basic equation of Fermi-Eyges theory. We recover their 
result that the transverse beam size at end-of-range is proportional to range, and their further result 
that (beam size)/(maximum size) versus (depth) /(maximum depth) is a universal function (which 
they express in closed form), independent of incident energy and stopping material. We generalize 
both results to arbitrary ion species, and we review the experimental data for protons and ions. 
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1 Introduction 



From time to time we teach a one week intensive course 'Techniques of Proton Radiotherapy' to 
persons involved in the design, upgrade or maintenance of proton radiotherapy facilities, or otherwise 
interested. An undergraduate degree in Physics or Engineering is required, but familiarity with 
particle physics is not. The course is equal parts basic physics, engineering (using the physics for 
useful calculations) and experimental devices and techniques. Lectures, a draft textbook, and some 
useful Fortran software and executables are available for free download [U [2] . 

These notes deal with deterministic transport theory, the topic (above all others) which time 
does not permit us to cover fully in the classroom. By this point in the course we have studied 
energy loss and multiple Coulomb scattering (MCS) of charged particles, and it is time to see what 
happens when those processes occur simultaneously with drift (motion along the beam axis). 




Figure 1: Top view of protons stopping in a water tank, with transverse scale exaggerated. Fermi- 
Eyges theory predicts the spatial and angular distributions of the full beam at any depth, and of 
the protons passing through the off- axis slit. 

FigureQ] shows a proton beam stopping in a water tank. The transverse scale is exaggerated 
(actual proton angles with respect to the beam are only a few degrees). The tracks were simulated by 
a simple Monte Carlo (MC) program. The first thirty were plotted without conditions. Thereafter, 
only those passing through a virtual off-axis slit were drawn. We'll sketch how the MC works, in 
order to contrast it with the main body of these notes. 

The MC divides the water into layers or slabs. Starting at the first slab, a proton is 'transported' 
as follows. Its incoming position and direction are projected to the midpoint, and its energy reduced 
by an amount appropriate to a half-slab, x and y deflections are chosen at random from a probability 
distribution that obeys the laws of MCS for water, that slab thickness, and that proton energy. The 
new directions are used to project the track to the downstream face of the slab, and another half-slab 
of energy is deducted. Given incoming position, direction and energy we now have outgoing position, 
direction and energy. We repeat for the next slab, and so on until the proton stops (energy reduced 
to zero). Then the whole thing is repeated for the next proton. 

Although real MCs such as Geant4 are more complicated, their basic principle is the same. Final 
quantities such as the rms spread of the beam at some depth, or the fLuence at some point, are 
found by scoring the MC results. Accuracy is limited by statistics, that is, by how many proton 
'histories' we have run. Thus the MC method is non-deterministic: its results change, consistent 
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with statistics, from run to run0 Unfortunately, Monte Carlo is the best way presently known to 
solve problems where the geometry through which the protons propagate is complicated. 

If, however, the geometry is simple (as in Figured]) a deterministic and much faster method 
is available: Fermi-Eyges (FE) theory, the main subject of these notes. Because of the restricted 
geometry, FE by itself solves relatively few practical problems. Among those are single scattering 
beam lines, the front end of double scattering beam lines, beam spreading in matter, and the 
omittance increase in a cyclotron energy degrader. Nevertheless it is worth studying because it is a 
building block in many other problems, such as pencil beam dose algorithms and finding the most 
likely path of a proton through a degrader. FE gives us a feeling for how a beam behaves in simple 
cases, before we tackle more complicated ones. 

2 Phase Space Diagrams 

As Figure [1] shows, relationships between proton positions and directions are complicated even in a 
simple situation. Phase space diagrams give us a systematic way of visualizing them geometrically, 
even if we have to bring in other methods to get quantitative results. In the phase space picture 
each proton propagates along the beam direction z and we examine its phase space variables 

x x' y y E 

at various 'depths' or values of z. x, y are transverse positions and x' , y' are slopes, which represent 
the proton's direction. In small angle approximation, which is always valid for radiotherapy protons, 
the slopes equal the projected angles 9 X , y . (When there is no possibility of confusion, we'll drop 
the subscript and refer to 9 X simply as 9. In other words, we'll use x', 9 X and more or less 
interchangeably.) E is the kinetic energyH A 'phase space diagram' is a scatter plot of x' vs. x, or 
y' vs. y, at a given value of z. For a sample, look ahead at Figure[3] Very soon, you will understand 
what is going on here and why such diagrams are useful. 

In the phase space picture E is a sort of 'ghost' variable, always present but not actually shown. 
It begins at the incident energy (say 160 Mev) and decreases as the protons progress down the 
beam line. When we start doing quantitative calculations we will need to keep track of it (that is, 
'transport' it) because — for instance — the strength of a given scatterer located at some z will depend 
on the energy of the protons at that z. In simple cases, E decreases as a function of z in more or less 
the same way for all protons, so each snapshot at a given z can be labeled with some E. In more 
complicated cases (for instance, the patient) protons with significantly different energies may reach 
the same point x, y, z because of heterogeneities and multiple Coulomb scattering. Such situations 
are beyond the scope of these notes. 

In some problems it may be more convenient to use a variable other than E as the 'longitudinal' 
variable. For instance, in problems involving magnets, the momentum pc (MeV) is convenient, 
whereas in dosimetry we may use the residual range in water (cm). Always, the longitudinal variable 
is some measure of the proton's speed. 

Phase space diagrams and the concept of the beam ellipse are much used in connection with 
particle accelerator and magnetic beam line design [5]. Their explicit use in papers on scattered 
beams is relatively recent, though even the early development of Fermi-Eyges theory implies a 
phase-space picture in the background. 

In oriented materials, scattering may be very different in x and y. We only consider homogeneous 
scatterers, where the scattering in x and y is the same, so it is sufficient to draw a diagram for one 
or the other. We'll use x, 9 (= 6 X = x'). That works until the point in the beam line where a beam 
limiting device (collimator) or something else breaks the symmetry between x and y. 

2.1 Model Beam Line 

Let's introduce phase space with a beam line of no practical use. It is a 'separated' beam line 
consisting of three thin scatters (MCS and energy loss only, no drift) separated by voids or 'drifts' 

1 To be sure, if the random number seed and everything else in the computer run is unchanged, we'll get exactly 
the same answer every time. Nevertheless, the statistical uncertainty is still there, but less obvious. 

2 In medical physics E is commonly used for kinetic energy, though T is sometimes used [3]. In general physics, E 
usually means total energy [4] as in E = mc 2 . We reserve T for scattering power. 
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Figure 2: Model beam line for phase space analysis: three scatterers separated by voids or drifts, a 
collimator with a central slit, and a final drift. 





Figure 3: Phase space evolution of the model beam (distributions at the entrance to each labeled 
element). For visual clarity the ellipse is drawn at 2.5 a rather than the conventional la; the 
standard ellipse would only contain about 390 of the 1000 phase space points. 
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(drift only, no MCS or energy loss), followed by an ideal collimator having a small slit on axis, 
followed by a final drift. In any real slab of matter MCS, energy loss and drift occur simultaneously. 
Eventually we'll take that into account. In the beginning it's easier to separate them. 

Figure^] shows proton rays in the model beam line with, as always, the transverse dimension 
greatly magnified. Figure|3] is a preview of the phase space diagrams corresponding to those rays. 
Each frame shows phase space going into that beam line element; thus S3 labels the phase space 
diagram obtaining at the entrance to the third scatterer. The next two sections explain these phase 
space diagrams. 

2.2 Effect of a Scatterer 

In a scatterer, points representing single protons receive random vertical kicks. The size of the kick 
has a Gaussian distribution with a width that depends on the scatterer material and thickness, and 
on the proton energy at that z. The distribution of kicks is independent of x for a uniform scatterer. 
The kicks may be positive, negative, small (more likely) or large. FigureS] illustrates all that. No 
point is untouched, except the lucky proton that does not scatter. However, because there is no 
drift in a thin scatterer, the horizontal position of each point is unaffected. 

2.3 Effect of a Drift 

In a drift, points in the upper half move to the right in proportion to their distance from the x 
axis whereas points in the lower half move to the left. In a pure drift (no scattering) the vertical 
position of each point is unaffected. The net result is a shear (not a rotation!) of the phase space 
distribution. Points on the x axis stay where they are. Figure[5] illustrates all this. To help we have 
drawn the proton trajectories (x vs. z) in the top part of the figure. 

2.4 The Beam Ellipse 

A subset of points in the phase space diagram can be surrounded by an ellipse. The entire behavior 
of the beam is encapsulated in the behavior of that ellipse, as we'll see. 

In scattering/drift beam lines (even when the two occur simultaneously) the ellipse has a rigorous 
quantitative meaning. It is an iso-density contour in phase space, as we'll show. That contrasts with 
accelerator theory (magnetic beam transport) where the beam ellipse is merely a convenience. It 
contains a defined fraction of the beam, but is not directly related to the phase space density. 

Figure[3] frame C (phase space at collimator entrance) shows a typical beam ellipse. The area 
enclosed by the beam ellipse is an important property of the beam known as its omittance. In 
the Appendix we show that the area enclosed by a tilted centered ellipse equals it x[ nt x max (or, by 
symmetry, n x' max Zint)H The area enclosed by a degenerate ellipse (point or line segment), and 
therefore the emittance of the beam so described, is zero. That means a beam can have a finite size 
(extent in x) or divergence (extent in x') or both, and still have zero emittance. Zero emittance 
means that position and angle are perfectly correlated, not that either is necessarily zero. 

2.5 Phase Space Diagrams for the Model Beam Line 

We are now ready to understand the phase space diagrams in Figure lBTl The first frame shows the 
ideal beam entering the first scatterer. The beam has neither spatial extent nor divergence, so all 
of the protons fall on a single point at the origin of phase space. (For better visual effect, we have 
drawn a small blob instead.) 

Leaving the scatterer into the first void, the beam has acquired a distribution of angles but still 
no size. Entering the second scatterer phase space has sheared. The beam now has both size and 
angular spread but x' and x are perfectly correlated and the phase space area (had the incident 
beam truly been a point) is still zero. 

3 That reduces to the more familiar (nab) if the ellipse is erect. 

4 In early dealings with phase space diagrams you will be tempted to think of them as cross sectional pictures such 
as might be obtained by putting a piece of film in the beam. That would in essence be a scatter plot of y vs. x. That 
is not what they are! They are (more abstract) plots of x' vs. x. 
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Figure 4: Phase space diagram of four protons acted on by a scatterer. 




Figure 5: Phase space diagram of four protons acted on by a drift. 
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Leaving the second scatterer, however, the phase space distribution will have an area no matter 
how perfect the incident beam. The exact correlation between angle and position is destroyed. 
'Angular confusion' has been introduced!! 

You can follow the rest by yourself: a shear, more scatter, another shear. We have drawn an 
ellipse around the beam entering the collimator|j 

Finally, a few of the protons pass through a narrow slit in the collimator. The rest stop. Protons 
emerging have no spatial spread left (except the slit width) but they have an angular spread. The 
rms value of that spread 9c is the quantitative definition of the angular confusion of the beam. Note 
that we have solved — by Monte Carlo simulation — a problem that might at first seem difficult: If 
an ideal incident beam makes its way down a sequence of scatterers and drifts, what is the angular 
spread of protons that happen to pass through a slit at the end? When we study Fermi- Eyges theory, 
we'll learn how to compute 9c another way, without simulation. 

Soon we'll also learn why we're interested in 9c ■ When we compute the transverse penumbra 
of a scattered beam (the sharpness of the shadow cast by a collimator and therefore, the sharpness 
of the dose distribution) 9c turns out to be the key and it, in turn, depends on the emittance as 
should be obvious from frame C of Figure[3j No emittance, no 9c- 

2.6 Emittance Change in a Drift 

Figure[7] shows a beam ellipse before and after a drift. x- m t does not change because points on the 
x axis do not move in a drift, a^ax does not change because there is no scattering. From our area 
formula we therefore conclude that the area enclosed by the beam ellipse, that is, the emittance, is 
conserved in a drift. This is a special case of Liouville's Theorem. 

2.7 Emittance Change in a Scatterer 

Consider Figure[3] panel VI. Scattering by itself does not increase emittance. If the incident beam 
were truly perfect, the phase space distribution in VI (degenerate ellipse) would still have no area. 
Emittance first comes about unavoidably in panel V2 just after the second scatterer. Two factors are 
needed for an omittance increase: the beam must have some size, or spread in x, and the scatterer 
must have some finite strength, that is, cause some random up or down kicks in 9. Indeed, we will 
show later that the increase (in quadrature) of the emittance in a thin scatterer equals exactly 7r 
times the product of the rms size of the beam at the scatterer times the rms angular deflection 
caused by the scatterer. A given scatterer will increase the emittance less, if it is further upstream 
where the beam is smaller. That becomes an important consideration in beamline design. 

2.8 Phase Space for a More Realistic Beam Line 

Figure[5] shows a modern double scattered passive beam lineQ We have drawn selected Monte Carlo 
generated rays, ones that pass through a virtual slit near each collimator edge. These justify our 
interest in 9c since the unsharpness of the collimator's shadow clearly depends on the rms cone angle 
of protons grazing the collimator. It also clearly depends on how far downstream the water tank is, 
showing why air gap is so important in proton radiotherapy. The dose edge is further degraded by 
scattering in the water tank, showing why depth in the patient is yet another important factor. 

Figure[5] shows the same beam in phase space, interleaved with projections onto the x axis to 
show fluence distributions in x. The first three panels are as before. The fluence entering the 
second scatterer is Gaussian. However, leaving S2 (entering the second air gap) two new things have 
happened. A collimator (not shown in Figure|8|) around S2 cuts off phase space at some radius, 
producing sharp edges. Also, S2 is stronger, and therefore the vertical kicks are greater, near the 
beam axis. The ensuing greater shear as protons drift to the patient collimator C is, by design, just 
enough to flatten the fluence distribution. 

5 The term 'angular confusion' was introduced by Andy Koehler in discussing protons [6], The fact that standard 
Fermi-Eyges notation for the same quantity is Qq [7] may be fortuitous. 

6 For visual effect we have drawn the ellipse much larger than would be usual, at 2.5 X a of the underlying Gaussian 
distributions, whereas the standard ellipse of which we will speak later would be at 1 X a. 

7 The first scatterer represents an upstream range modulator stopped at one of its M thickness steps. To solve the 
full problem, we'll eventually need to solve it for each step separately, then sum over M using appropriate weights. 
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Figure 8: Double scattered beam line with selected trajectories of protons that graze the downstream 
edges of the patient collimator. SI represents one step of a range modulator. 
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Figure 9: Phase space diagrams for the double scattered beam with x projections interleaved, (n.b. 
The setup used for these differs in two points from Figure[8l There is a collimator surrounding S2, 
and the 'patient collimator' is a half-beam block with one edge on axis rather than the symmetric 
one of Figure[5]) 
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The patient collimator chops off protons on the left and right creating sharp edges at A3. Those 
edges shear in the ensuing drift so their projections entering the water tank smear out, causing 
penumbra which depends not only on the angular confusion (vertical extent of the phase space 
distribution) at C but also on the length of the drift (air gap). Multiple scattering and drift in the 
water tank blur the edges still further. 

2.9 Summary 

The phase space coordinates of a proton at depth z are x, x' , y, y' and E. 

A phase space diagram is a scatter plot of x' vs, x. Each point represents one proton. 
The diagram is not a cross section of the beam! 

In a thin scatterer each point acquires a random ± vertical kick, whose rms size depends 
on the strength of the scatterer. Its x position is unchanged. 

In a pure drift, the phase space distribution shears, points above the x axis moving right 
in proportion to their distance from the axis, points below moving left. Points on the x 
axis do not move. 

A phase space distribution can be characterized by an ellipse which is an isodensity 
contour (to be shown). 

The beam emittance is the area enclosed by the ellipse and equals 7r x' int x max = ir x' max x; nt 

Emittance is conserved in a drift, and increases (in quadrature) by ir x rms 8 n in a thin 
scatterer (to be shown). Thus a given scatterer has greater effect, the larger the incident 
beam. 

The evolution of the beam ellipse (three parameters) can be computed by Fcrmi-Eyges 
theory (to be shown). That amounts to complete knowledge of the beam at every z 
provided the requirements of FE theory are met. 

3 Miscellaneous Topics 

This section gets a few issues out of the way before we begin Fermi-Eyges theory proper. 
3.1 Review of Gaussians 

The standard ID Gaussian probability density P and its normalization, mean and variance are 

1 ( U — Uq \ 2 



1 / u - Up \ - 
2\ a ) 



P(u) = -=^e 2 V a ) (1) 
V27T a 



/oo 
P{u)du = 1 (2) 
-OO 

uP{u)du = < U > = Uq (3) 

/oo 
u 2 P{u)du = <u 2 > = cr 2 (4) 
-oo 



/. 



P(u) stands for a probability density whose dimensionality is indicated by the number of arguments. 
Thus P{x,9)dxd9 is a joint probability, P{x\9)dx is a conditional probability and NP(x,8) is a 
number density in phase space given N incident protons. 
The 2D 'cylindrical' Gaussian probability density is 

1 / r \ 2 

P(r,<f>) - -^e~ (5) 
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which is also normalized 

i>2tt />oo 



JO 

However 

/■27T />0O 

and Eq.[5] is often written 



o Jo 



P{r,4>)rdrd(j} = 1 (6) 
r 2 P{r,(j))rdrd(j) = 2r^ (7) 



P(r,0) = — ^ e (8) 

When reading a paper, consider whether the authors use Eq.[5]or Eq. [5] because of the y/2 difference 
in the Gaussian width parameter (ro or oy). 



3.2 The Gaussian Approximation to MCS 

Angular and fluence distributions in Fermi-Eyges theory are Gaussian: FE alone can never predict 
a non-Gaussian distribution of anything. That flows from the basic assumptions (infinite uniform 
slabs) and the fact that MCS is very nearly Gaussian. How good is that approximation? 

FigurcfTOl shows the outgoing projected angle distribution of 158.6 MeV protons traversing 1 cm of 
water, according to three variants of MCS theory^ 'Highland' is a Gaussian with a width parameter 
#o given by Highland's formula, a simple and quite accurate parametrization of the full theory. 
'Hanson' is also a Gaussian, but its width parameter is obtained by computing the full theory, then 
finding the Gaussian that best approximates it. Finally, 'Moliere' is the full theory of the angular 
distribution. It has been compared to experiment and is known to be correct to a few percent [8]. 

Figure [10] shows that the central part of the angular distribution is indeed Gaussian out to 
about 2.5cr, and that the three methods agree very well in that region. The ±2.5 a region contains 
about 99% of the protons according to the normal probability integral [9]. If the projected angle 
distribution is Gaussian, then in small angle approximation the projected spatial distribution on a 
screen or measuring plane somewhere downstream will also be Gaussian. That justifies a theory 
that predicts nothing but Gaussians in the very simple cases for which it holds. 

Most problems in proton radiotherapy are treated adequately in the Gaussian approximation, 
with a width parameter given by one or another version of Moliere theory. In the rare case where 
that is inadequate, the problem is usually not the 'single scattering tail' of the Moliere distribution 
evident in Figure [TD], but the halo of charged secondaries from non-elastic nuclear reactions that 
soon appears around any proton beam traversing matter [TO] . 



3.3 Relativistic Single Particle Kinematics 

Any introduction to special relativity will provide the basis for the following equations. If E is the 
kinetic energy of a particle we define its reduced kinetic energy as 

t = E/mc 2 (9) 

where mc 2 is the particle's rest energy. Then j3 = v/c is given by 

? = (10) 

where v is the particle's speed and c is the speed of light. The kinematic quantity pv, common in 
MCS theory, is given by 

pv = T -±1e (11) 
r + 1 

where p is the momentum. For completeness, the momentum itself is given by 

{pcf = (T + 2)mc 2 E (12) 

which we would need were we dealing with charged particle transport in magnets. These three 
equations are convenient because their relativistic and nonrelativistic limits are obvious at a glance 
and because they avoid small differences between large quantities which can occur in relativistic 
kinematics if one is not careful. 



For details, please see [8]. 
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Figure 10: Projected angle distributions for 158.6 MeV protons traversing 1 cm of water. On a 
linear plot the Moliere/Fano distribution is indistinguishable from Gaussians using the Hanson or 
Highland 6q. On a log plot, the correct distribution peels away at 2.5a, and is more than lOOx 
higher at 5a. 

3.4 Completing the Square 

The technique called 'completing the square' allows us to evaluate certain definite integrals. The 
formula is 

au 2 + bu + c = a(u - h) 2 + k (13) 

where 

b b 2 
h=-7r , k = c~— 14 
2a ' 4a ' 



3.5 Scattering Power 

An important concept in transport theory is the scattering power 

m = d5) 

dz 

T, the rate of change with z of the variance of the projected MCS angle, bears a superficial resem- 
blance to stopping power S = —dE/dz, the rate of energy loss. However, the analogy is imperfect. 

S, for which there is one generally accepted theory [3], is a function of the proton's speed or 
energy and of atomic properties of the stopping material at the point of interest. If we compute 
S everywhere in a finite slab and integrate, we obtain the correct energy loss in the slab. But if 
(following Rossi |11| ) we compute T as a similarly 'local' function of speed and stopping material 
and integrate that, the answer is substantially too large, particularly for thin slabs. 

That might not seem to be a problem because the generally accepted theory of multiple Coulomb 
scattering — Moliere theory [T2J[T3J[T3] — gives 9 X accurately for arbitrary ions, materials, slab thick- 
nesses, and incident energy. But Moliere theory deals only with finite slabs and is therefore not 
directly suited to transport calculations, deterministic or Monte Carlo, which require a formula for 
the instantaneous rate of change of 9 X . T is most usefully understood as satisfying that requirement: 
a differential approximation to Moliere theory, or a function which, when integrated, reproduces the 
Moliere/Fano/Hanson angle [8] more or less accurately. 



12 



We find that the formula for T must be nonlocal. The rate of change of X with respect to z for 
(say) 50 MeV protons in water depends to some extent on how much water they have already gone 
through. There are different ways of expressing that nonlocality [TB] . 

We repeat here, from [T3], a few formulas necessary later. Our own 'differential Moliere' scattering 
power is 

E S V 1 



TdM = f&M{pw\,pv) X — (16) 

V pv J X s 

where E s — 15.0 Mev. The correction factor, which measures nonlocality by the decrease in pv, is 

f dM = 0.5244 + 0.1975 log 10 (l - {pv/p^f) + 0.2320 log 10 (pu/MeV) 

- 0.0098 log 10 (p V /MeV)log 10 (l - (pv/ PlVl ) 2 ) (17) 

PiVi is the initial value and pv is the value at the point of interest. The scattering length Xs (cm) 
is a property of the material at the point of interest given by 



1 „ 2 z 



aNr 



{21og(33219(AZr 1/3 ) - l} (18) 



where p is density, a is the fine structure constant, N is Avogadro's number, r e is the classical 
electron radius and A, Z are the atomic weight and atomic number of the scattering element. In 
compounds or mixtures, atoms act independently ('Bragg rule') and the slab is equivalent to very 
thin sheets of each constituent in the correct proportion. That picture leads directly to 

^-?-(^), (19) 

where Wi is the fraction by weight of the i th constituent. Xs is a material property very similar to 
Xq, the radiation length. Table [1] compares Xs with Xq for some materials. (Note that it is mass 
scattering length and radiation length that are tabulated.) 

When integrated, T<jm reproduces Moliere/Fano/Hanson theory, as well as experimental data at 
158.6 MeV [8 , to a few percent over a wide range of materials for normalized slab thicknesses Az/R 
from 0.001 to 1 [15]. There is indirect evidence that it also works for mixed slabs [15], but that has 
not yet been tested experimentally for T^m or any other scattering power. 

The nonlocal correction to T seems paradoxical, implying that the interactions of a proton in 
water at 50 MeV depend upon how much water has been traversed. However, multiple Coulomb 
scattering is not a primitive process! It can be viewed as a statement about the statistics of a large 
cohort of protons each of which has suffered a large number of atomic encounters. Without the 
nonlocal correction factor their rms spread in angle would be exactly proportional to the square 
root of slab thickness, but it is not, because MCS is not exactly Gaussian. Therefore the degree of 
approach to 'Gaussianity' also plays a role. That is a statistical statement with no implications for 
the underlying single scattering process which indeed depends only on a proton's speed and impact 
parameter at the point of interaction, and on the atomic properties of the scattering material. 



4 Fermi-Eyges Theory 
4.1 History 

Fermi (concerned with the propagation of cosmic rays) considered the joint probability of position 
and angle, due to MCS, of a single charged particle entering a uniform medium (the atmosphere). 
He evidently did not consider his derivation worth publishing, since Rossi and Greisen [16] state in 
a footnote 

'The developments in this article follow closely a lecture given by Professor Fermi at the 
University of Chicago in the summer of 1940 and include some unpublished results. The 
writers wish to express their sincere appreciation to Professor Fermi for allowing them 
to make use of these results.' 
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Later, Eyges [17] included the effect of energy loss, which Fermi had ignored. Still later, the theory 
was generalized to cover a non-ideal (but still Gaussian) incident beam, a stack of slabs of different 
materials, and improved formulas for the scattering power. We can regard this generalized theory 
as a way of propagating the beam ellipse (hence the phase space distribution) through homogeneous 
mixed slabs, and we'll refer to it simply as Fermi-Eyges or FE theory, dropping the 'generalized'. 

Here are a few comments for the student reading the early papers to reconcile those with our 
version. Neither Rossi and Greisen [16] nor Eyges [17] use the term 'scattering power'. However, 
Rossi's book [11] uses 6 2 ~rad 2 /(g/cm 2 ) and clearly recognizes that some formulas for 9 2 are better 
than others, as well as the issues raised by more accurate theories such as Moliere's. The scattering 
power more or less built into Eyges' derivation (our Tfr) is the worst possible choice [18j [T5] . The 
term 'mass scattering power' was introduced by Brahme |19] who was also an early user of FE 
theory in radiotherapy (electrons). The same paper gives a way (not generally used nowadays) of 
incorporating Moliere theory into the scattering power. 

Early papers alternate between space and projected angle. We use projected angle exclusively, 
which makes our scattering power half as large. Also in early papers, depth is commonly expressed 
in units of the radiation length of the material, which is awkward if there is more than one material. 
For the same reason we normally use scattering power T ~ (rad 2 /cm) rather than mass scattering 
power T I p ~ (rad 2 /(g/cm 2 ), range R rather than mass range pR, and so on. 



M l -2 
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i; 
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Ay. 



Figure 11: Stack of three homogeneous slabs for Fermi-Eyges discussion. Pb/Lexan/air simulates 
the upstream end of a double scattering beam line. 



4.2 The Basic Theory 

Consider Figure QT] Our eventual goal is to find the phase space distribution and related quantities 
at any depth z if an arbitrary Gaussian beam enters at z = 0. Let z be the beam direction and 
x, y transverse directions (y is up and x,y,z form a right handed frame). For the present, suppose 
a single proton enters the first slab, whose upstream face is at z — 0, along the z axis. In a very 
short but mathematically sophisticated paper, Eyges [17] showed that the probability of finding the 
proton later at some z > with x in dx and 9 in dO is 

1 A x 2 - 1A x xQ + A 2 2 

P(x,0)dxd6 = — l —e~2 B dxd6 (20) 

27rv B 
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A Q {z) = 


/ T(u) du 




Jo 


Ariz) = 


/ {z — u) T(u) du 




Jo 


A 2 (z) = 


/ (z-u) 2 T{u)du 




Jo 


B(z) = 


A A 2 - A\ 



The A n are moments of T namely 

(21) 
(22) 
(23) 

Jo 

and 

(24) 

From now on, we suppress the z dependence, it being understood that the A n and any related 
quantities are functions of depth. 

B is greater than zero if there is any scattering at all, as the following proof by Jette [20 shows. 
The first line is true by inspection for any z > because T > (the mean squared angle can only 
increase with z). 

< A [ T(u)((z-u)-^-) 2 du 
Jo ^ Aq/ 

= Aq f Z (t{z - uf - 2T(z - u)4^ + T(^) 2 ) du 

JO V Aq Aq / 

= A (A 2 -2A 1 ^ + A (^f) = A 2 A -A 2 = B 



To find the distribution of 9 irrespective of x we integrate P(x, 9) over x. Completing the square 
in Eq. (|20p by using Eq. (|13[) with u = x we find 

P( X( 0) = —L= e 2Aq e 2B (X Aq (25 ) 
27rv B 

l ° 2 iA q A 1 2 

/OO i />OC ($ — \j\ 
P(x, 9)dx = = e 2A o e 2 B a q dx 
2nVB ' 



— OO 



The term Ai9/A , independent of x, is equivalent to a shift in origin which we can ignore since the 
integral runs from — oo to oo. Using Eq. we find 

P{9) = -y^=e 2A (26) 

Comparing with Eq. (QJwe find that Aq equals the variance of 9 

Aq =<9 2 > (27) 

which was already obvious from Eq.[2l] and the definition of T. 

By virtue of the symmetry of Eq. (1201) any equation in 9 can be turned into a corresponding 
equation in x by swapping x, 9 and Aq,A 2 . Thus Eq. ([23)1 becomes 

P(x y Q) = —L= e 21 2 e 2B ( A 2 (28 ) 
27rv B 

while Eq. (l26l) becomes 

_lx^_ 

PW = ism r ^ (29) 

and A 2 equals the variance of x. 

A 2 = < x 2 > (30) 
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That is less obvious, but a simple physical argument will be given later fSection l5.1|) . 

In addition to the gross variance of x and 9 we require 9 C , the variance of 9 at a given x (the 

angular spread of protons emerging from a narrow slit at some x) and x 2 s , the variance of a; at a 
given 9 (harder to visualize) According to probability theory the conditional probability of 9 given 
x is 

P{ e\x)d8 = P ^\ dxdd (31) 

v 1 ' P(x)dx v ; 



2ttB 

showing that, given x, the most probable 9 is 



IA 2 A 1 2 

— r 2 d6 (32) 



2 



Furthermore 



9 2 



/oo r> 
tf 2 P(0|z)^ = — (34) 
-oo A 2 

9c is independent of x: the angular spread of protons emerging from a slit is independent of the 
tranverse position of the slit.F"! 



Similarly 



>2 



P(x\9)dx = \I^b e 2 B A ° ' dx (35) 



showing that, given 9, the most probable x is 



Ai 



(36) 



and 



B 



xi s =<x\i>= x<P{x\6)dx = — (37) 

which is independent of 9. The transverse spread of protons at a given inclination to the axis is 
independent of the inclination. 

We have given simple physical interpretations of A$ and A2 . To interpret A\ consider < x 9 > , 
the covariance of x and 9. The average can be arranged in either of two ways. Choosing the outer 
average to be over x (denoted <> x ) and denoting l 9 given x' by 9\ x we have 

< x 9 > x = < x x 9\ x > x 

= <x<9\ x > g > x 

= < x((A 1 /A 2 )x) > x 

= (A x /A 2 ) <x 2 > x = (A 1 /A 2 )A 2 = Ai 

The third line follows from Eq. (f33|) and the fact that the distribution of 9\ x is Gaussian so the mean 
equals the most probable value. The fourth line follows from Eq. (J30j) . Thus 

Ai =<x9> (38) 

Of course, the same result is obtained doing the averages the other way. 

So far we have linked Aq to the gross angular spread (Eq. [27| . A\ to the correlation between x 
and 9 (Eq.l38|) and A 2 to the gross transverse spread (Eq.l30j). Further interpretations of the A n and 
B, of a more geometric sort, will appear in the next section. As an application of the ones given so 
far, consider 

c A 2 A 2 A 2 < x 2 > y ' 



1)2 



9 The name x e ff is explained later. 

10 That is somewhat obvious from first principles. Although the mean angle of protons emerging from a slit depends 
on its transverse position, the 'scattering history' of those protons, in small angle approximation, is statistically the 
same as on-axis protons and therefore the spread about the mean is the same. 
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The angular spread at a point (spread of protons emerging from a narrow slit) is less than the gross 
angular spread unless < x 9 > = 0. 

Note: many writers absorb the physical interpretation of the A n into the notation, writing 
Eq. (|20p for instance as 

1 Wx 2 - 2x~3x6 + ~x 2 9 2 

P{x,0)dxd0 = 1 2 x 2 9 2 -^f dxd9 

2n\Jl^lP-xf 

We prefer Eyges' notation as being easier to read. It also emphasizes that the A n are simply three 
numbers that can be calculated as a function of z, the depth in the stack. 



4.3 The Beam Ellipse 

The contour in phase space on which the probability density Eq. (|20|) equals e -1 / 2 = 0.606 = 61% 
of its central value is found by setting 

F{x,6) = A x 2 - 2Aix9 + A 2 6 2 = B (40) 

According to the analytic geometry of conic sections Eq. P0|) describes an ellipse because B > 0. It 
is a centered ellipse because F is invariant under x — > —x, 9 — > —9 (no linear terms in x or 9). 

The extrema 9 of the ellipse are found as usual by setting dF/dx equal to 0, giving x — Ai9/Aq. 
Putting that into Eq. (|40|) and simplifying we find that 



9 = ±y/A, =<9 2 > 1/2 occurring at x = ±A 1 /^f~A^ (41) 

Similarly 

x = ±\f~A 2 = < x 2 > 1/2 occurring at 9 = ±Ai/\/A^ (42) 

In other words, the rms values of x and 9 define the bounding box of the ellipse. Figure[l2] shows 
what we have learned so far in A, B notation. Figure H31 shows the same thing in x, 9 notation. 
In Appendix [Al we show that the area bounded by a centered ellipse in x,y is 

7T ?/int ^max — 7T ^int 

According to Eq. (1341) 9c is the rms spread in 9 at any given value of x, including 0, which means 9c 
is the y intercept of the ellipse. Eq. (|23|) gives the conjugate maximum of the ellipse, x ma , K — V 'A 2 . 
Thus _ 

e = 9 C <x 2 > 1/2 = ^B/A 2 ^/T 2 = \fB (43) 

or 

area bounded by the ellipse = it Vb (44) 

which is the interpretation of B. In accelerator physics the bounded area is called the emittance of 
the beam [S]I]3 However, in accelerator physics the beam ellipse is merely a convenience (we know 
how to transport it through magnets and drifts) and Gaussian phase space density is usually only 
a rough approximation. By contrast, in proton transport through matter the phase space density 
really is Gaussian to a very good approximation and iso-density contours really are elliptical. 

A\ = < x 9 > is related to the tilt of the beam ellipse. It is conventional [7] to derive an equation 
for VP, the angle from the x axis to the principal axis of the ellipse, namely 

2 A 

tan 2* = — (45) 

-42-^0 

This must be used with discretion as the denominator warns us: A 2 and Aq have different dimensions 
and cannot be added. '5 refers to the angle as measured in a particular drawing of the beam ellipse, 
and the A n in Eq. (|45|) are appropriately scaled versions of the real A n . 

From Eq. (|3"2"j) the most probable 9 given x is 9 = (Ai/A 2 ) x, a line which we can identify with 
OR in FigurcQ21 Similarly from Eq. (|55|) we can identify OQ as the locus of the most probable x 
given 9. 



11 Confusingly, in the early accelerator literature area/V (our yf~B) is sometimes called the emittance. If we quote 
the emittance as e.g. '10 tt mmmrad' it is clear we mean area, with y/B = 10 mm mrad. 
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Figure 12: The beam ellipse in A ni B notation. 
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area enclosed = n (<6 2 ><x 2 > - <x0 > 2 ) y = 

Figure 13: The beam ellipse in x, 9 notation. 



4.4 Drawing the Ellipse Given the Moments 



Instead of deriving ^ let us solve a related problem: Draw the ellipse corresponding to given Aq, A\ 
and Ai. In a computer program, any ellipse is most conveniently drawn using the parametric form 
of Appendix [A] whose equations Eq. (|87)l and Eq. (|88|) we can rewrite 

x = x cost (46) 
9 = 9 sin(< - 5) (47) 

x and 9 define the bounding box in graph units. To generate points on the ellipse, let t range in 
steps from to 2-k. To find <5 note that Eq. (|102j) . in our present notation, yields 

cos5 = e = = V^o^ - A\ 



6 = siiT 1 (A 1 /^A A 2 ) (48) 

Unlike using B (always positive) as the third ellipse parameter, this procedure preserves the sign 
information (converging or diverging beam) in A\ because the principal range of sin -1 is —n/2 to 
7r/2. To plot the ellipse in (mm, mrad) if y/A~2, \/Ao are in (cm, rad) use 

x = 10 V^b 
6 = 1000 

in Eqs. (ggj and fljZj) , 

It may be (cf. FigureE]) that we do not know the A n explicitly but wish to draw an ellipse around 
a set of points obtained by Monte Carlo, or otherwise. In that case we find the statistical rms values 
of Xi and 9i and fit a line to 9i(xi) to find < d9/dx >. Then 

a = tan- 1 (< 9 2 > 1 ' 2 / < x 2 > 1 ' 2 ) 

tp = tan _1 (< d9/dx >) 

S = sin _1 (tan(2'i/;)/tan(2a)) 

In Figure[3]we multiplied < x 2 > 1 / 2 and < 9 2 > x / 2 by 2.5 to obtain the larger ellipse. 



4.5 Ellipse Examples 

FigureQT] is a reality check which puts to work everything we have learned so far. It shows beam 
ellipses corresponding to a 127 MeV proton pencil beam propagating in a near stopping thickness 
(0.97 x range) slab of water, divided into five equal parts. We used T^m for the scattering power. The 
bounding box of the final ellipse is shown. Its x rms agrees well with the measurement by Preston 
and Koehler [2"T] . 

Next, for a problem of practical importance that can be fully handled by Fermi-Eyges theory, 
let's return to Figure [TT1 which represents one step of a compensated range modulator [I] and the air 
between it and the second scatterer. Anticipating the next section, Figure[T5] shows the evolution 
of an ideal 230 Mev incident proton beam in Pb/Lexan/air. The Pb approximates a thin scatterer. 
Only the distribution of 9 changes. The Lexan introduces more scattering as well as some drift 
making the emittance non-zero. The air (here divided into three slabs) is an almost pure drift: no 
increase in 9 rms or emittance. Now let's see how we compute those ellipses. 



4.6 Transporting the Beam Ellipse Through a Slab 

We want to generalize from a perfect beam entering a single slab to a beam of finite emittance 
entering a stack of slabs (each homogeneous and transversely infinite) . Some slabs may be air (little 
scattering) or void (no scattering). We call this 'mixed slab' or 'stack' geometry. 

Consider a perfect beam entering a homogeneous slab extending from to zi with an imaginary 
break at z\. Take (for instance) the second moment at z 2 namely 

f 2 \zi- z') 2 T{z')dz' = f 1 \z 2 - z') 2 T(z')dz' + f 2 \z 2 - z') 2 T(z')dz' 

JO JO J zi 
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Figure 14: Evolution of the beam ellipse of an ideal 127 MeV proton beam in near stopping thickness 
(0.97xrange) water divided into five equal slabs. The dashed lines represent the uncertainty range 
of a measurement by Preston and Koehler |21) . 
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Figure 15: Evolution of the beam ellipse for a 230 MeV proton pencil beam entering a Pb/Lexan/air 
(respectively 0.234/14. 6/85. lcm) stack. The air is divided into 3 equal slabs. 
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Consider the first integral on the RHS. If we add and subtract z\ we may write it 



Z \(z 2 -z 1 ) + (z 1 -z')) 2 T(z')dz' 







zi 



{{z 2 - Zl ) 2 + 2(z 2 - zi)( Zl - z') + (zi - z') 2 ) T(z') dz' 
= {z 2 - zx) 2 A 0l + 2(z 2 -z 1 )A n + A 21 
where Aqi stands for the zeroth moment at z\ and so on. Thus 

\z 2 - z') 2 T{z')dz' = (z 2 - Zl ) 2 A 01 + 2{z 2 - Zl )A n + A 21 + f * \z 2 - z') 2 T(z') dz' 



Since z\ and z 2 are perfectly general we can set z\ to and z 2 to z. Proceeding similarly for the 
other two moments and letting unprimed denote a moment at and primed a moment at z we find 



A' Q = A + / T(z')dz' (49) 
Jo 

A[ = A x + A z + [ (z- z') T(z') dz' (50) 
Jo 

A' 2 = A 2 + 2A lZ + A z 2 + ( (z- z') 2 T{z')dz' (51) 



There is no physics in Eqs. (14"9T l5"T j) . They express a mathematical property of moments based on 
the additivity of definite integrals. The initial moments need not come from scattering. They may 
represent the parameters of an incident beam. Any set is valid provide Aq > 0, A 2 > and 

B = AqA 2 - A\ > 

If Ax — < x 9 > is negative the beam ellipse slopes down (9 decreases as x increases) which 
(recalling our discussion of phase space diagrams) describes a converging beam. That cannot be 
produced with scatterers, but can be, and often is, produced with magnets. After a sufficiently long 
drift, any converging beam will diverge as the protons cross the beam axis. That too is built into 
Eqs. (|39>[5T). 

z' stands for position along the stack. Its role in expressions like (z — z') 2 is obvious. In T(z') it 
is more subtle. It is a bookeeping tool meaning 'where we are in the stack'. As shown by Eq. (|16l) . 
one possible choice for T, we will need a kinematic quantity pv (which depends on kinetic energy, 
which depends on z') and atomic properties of the current material (which also depend on z'). The 
beam transport program will have to include a procedure which accepts z' and, knowing the incident 
energy and how the stack is constructed, returns the variables upon which T directly depends F^l 

To summarize this section, we have derived three equations which transport an initial set of A s 
(Fermi-Eyges moments or ellipse parameters) through a finite homogeneous slab. The initial ^4s 
can represent either an incident beam or previous scattering. If they are zero the incident beam 
is ideal and we recover the original definitions of the moments. The equations can be repeated as 
often as necessary to transport the beam through a stack. The kinetic energy, or some equivalent 
longitudinal quantity, must also be transported, of course. 



4.7 Emittance Change in a Drift 

If a slab is void T(z') = 0, the integrals in Eqs. (|49l) through (I5T1) vanish, and the remaining terms 
transport the incident ellipse through a drift. By direct calculation 

B' = A' A' 2 - A\ 2 = A A 2 - Ai 2 = B 

for any z, recovering the fact that emittance is conserved in a drift (despite the fact that Ax, A 2l 
the phase space distribution, and the beam cross section all change, except in trivial cases). 

12 In our implementation that subroutine, in module SPWRsubs.FOR, is called WhatsHere(,z', . . .) . 
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4.8 Emittance Change in a Scatterer 

Eqs. (H§l - l5Tj) have another important consequence. Suppose a given beam (given As) falls on a thin 
scatterer so that we can assume T is nearly constant and evaluate the integrals accordingly. Then 

A' a = A + TAz 

A[ = A-l + A Az + f (Az) 2 /2 

A' 2 = A 2 + 2Ai Az + A (Az) 2 + f (Az) 3 /3 

where T is the value at the midpoint of Az. To lowest order in Az 

B' = AqA 2 — A'i ~ A A 2 - A\ + TAzA 2 

or 

AB = A{( 2 ) = A 2 (TAz) = <x 2 >xA<9 2 > 

If we take the square root, the quadratic difference on the RHS is just 9 , the rms strength of the 
thin scatterer [5]. Thus the quadratic change in emittance in a thin scatterer is 

■k\J e| - e\ = 7r x rms 9 (52) 

where x rms is the rms size of the beam impinging on the scatterer. If a beam of known emittance 
enters N scatterers separated by drifts, the final emittance is 7re where 

N 

g2 = e bcam + x rms,i X @0,i 

Thus, if a passive beam line can be approximated by thin scatterers separated by drifts (as is 
usually the case) the final emittance can be approximated without the full Fermi-Eyges treatment 
by just computing the beam size at each scatterer, which is easier (see Section lBTTj) . Of course, each 
9o must be computed with due regard for the proton energy at that scatterer. 

Eq. (|53[) has important consequences for beamline design. As Figures[5]and[S]show, the transverse 
penumbra of a double scattered beam depends on the angular confusion at the patient collimator, 
9c, which is proportional to e (Eq.l43p. Therefore the penumbra can be reduced by placing the 
second scatterer further upstream, where the beam is smaller. (That is also obvious geometrically 
from Figure|8l) 

If the incident beam is small, the first scatterer (no matter how strong) will cause relatively 
little emittance increase. That accounts for the small penumbra of single scattered beams. In fact, 
their penumbra is usually dominated by the air (which acts as a second scatterer) between the first 
scatterer and the patient. 

4.9 Differential Form of the Transport Equations 

Broadly speaking there are two applications of transport theory: beamline design and dose recon- 
struction in the patient. Passive beamlines are composed of slabs (though some may be nonuniform) 
so Eqs. (|49l - l5"Tj) are a good starting point. In the patient, however, material properties change voxel 
by voxel and a finer-grained approach is needed. We wish to compute the changes in moments as the 
beam passes through a voxel of length Az, assuming the material is known. These are easily derived 
from Eqs. (|49l - l5"T]) by assuming T approximately equals its central value T in each voxel ('midpoint 
rule'). If, following Kanematsu [18J . we rewrite the result to minimize the number of multiplications 
we find 

AA = TAz (54) 
AAi = (A Q + (f/2)Az)Az (55) 
AA 2 = (2Ai + (Aq + (f/3)Az)Az)Az (56) 

These equations are usually embedded in a general procedure which also transports other quantities 
(e.g. kinetic energy or pv or residual range) through Az. 
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4.10 Equivalent Sources 



It is sometimes convenient to mentally replace the beamline upstream of some measuring plane zmp 
by an equivalent source. For example, we might want to replace the Pb/Lexan/air stack which 
represents one modulator step of a double scattering system by an equivalent point source some 
distance S upstream of the location of the second scatterer. Knowing S gives us the effective 'throw' 
of the first scatterer, allowing us to adjust its strength to produce the desired Gaussian width on 
the second scatterer. 

There are three different ways of defining an equivalent source. The first two are in ICRU Report 35 
[7]. Our formulas agree with that reference, of course, but our derivations are more geometric. 



4.10.1 Effective Extended Source 

The effective extended source is that ellipse which drifts into the given ellipse at zmp in a drift 
distance S c s- Since there are infinitely many such ellipses, we require in addition that the source 
ellipse be erect, < x9 > = A\ =0. FigurefTBl (top) illustrates the situation. Consider point Q. The 
drift distance of a general phase space point (x, 9) from an initial (0, 9) is 

S = \ (57) 
From the known phase space coordinates of Q fFigure [T2"]) we therefore find 

A 1 j^M A 1 A x 

S off = -=— = — whence z cS = z M p - ~r (58) 

VA, A A 

The effective extended source reproduces the given ellipse exactly. Emittance is the same since 
omittance does not change in a drift. Since points on the x axis do not move in a drift, the rms size 
of the (erect) source ellipse is 

x eS = JWjA (59) 
as shown earlier fEq.[37|). That justifies the notation x D g. 



4.10.2 Virtual Point Source 

The virtual point source is the ellipse that would drift into the line OR, the locus of the most probable 
proton direction given x. Therefore it is the point from which the protons observed at zmp seem 
to emanate. For a unique definition we again require that the source ellipse be erect (Figure[T6l 
bottom). A line is a degenerate ellipse (emittance = s/B = 0) and an erect line is a point source 
with finite angular spread. Its parameters can be read from Figure flTjl Using Eq. (|57[) and the 
coordinates of R we find 

^ A2 v, A2 fnn\ 

whence z v ir = z mp ~ -r~ (60) 



v " A./JM A, ~ v " — A, 

The gross angular spread of the given ellipse is the quadratic sum of the spread of the virtual point 
source and the angular confusion: 

¥ = 9^ + el (6i) 

as can be verified by expressing each term in terms of the A's. 



4.10.3 Effective Scattering Point 

The author, ignorant at that time of Fermi-Eyges theory, defined [5] a third equivalent source 
distance using the back projected asymptote of u x {z) fFigure lTTl) . Kanematsu later called this the 
effective scattering point [22) . It can be identified with the upper right corner of the bounding box. 
Using Eq. (T57| as before we obtain 



Scsp = -= = \l — whence z csp = z M p - \j — (62) 
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Figure 16: Ellipses for the (top) effective extended and (bottom) virtual point sources. 
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Figure 18: Source locations v. depth of measuring plane using a 12 cm Lexan slab starting at z = 0, 
followed by air, and an ideal 160 MeV proton beam. Full circles use T^m ', dashed lines use Tic [E] ! 
open circles (zbg) use Highland's formula [8 (compare with z csp , full circles of same size). 




Figure 19: Practical application of the effective scattering point. The first scatterer of the Burr 
Center 'STAR' beam uses Pb and Lexan degraders controlled by a computer to produce the desired 
SOBP. As the degraders change so does the effective scattering point (filled squares). Note the jump 
at the major binary transition. The related change in 1/r 2 is compensated by adjusting the dwell 
time of each configuration in the beam. 
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Looking at Figurc[l2] and keeping Eq. (|57|) in mind we can immediately write 

SVir > S csp > Ses whence z vir < z csp < z cS (63) 

After any substantial drift B, though it remains constant, is much smaller than its largest possible 
value AqA 2 (see Figure fTS"]) . The three source distances become very nearly equal, and the potentially 
nettlesome question of which one to use in a given problem becomes moot. 



4.10.4 Examples 

Suppose we have a system not unlike the Burr Center eye treatment line, a 160 MeV proton beam 
traversing 12 cm of Lexan followed by sa 90 cm of air. (The residual range is 3.9 cm H2O.) FigurefTgl 
shows various computations of the equivalent source as a function of the depth of the measuring 
plane (MP). Full lines, full circles are computed with the preferred scattering power T^m', dashed 
lines are computed with Tic j the same as T^m without the single scattering correction [15] . The 
difference is minor. Open circles represent the method described in [8] (back projection of the beam 
envelope at the la level) which we often use in place of Fermi-Eyges theory. In principle it should 
be compared with z esp , but it uses Highland's formula, making the comparison questionable. 

Air has a noticeable effect on the calculation. The emittance increase in air is s» 2.4x. If 
we substitute helium or hydrogen, the emittance increase is much less and the curves rise almost 
vertically (little change in z xxx with zmp)- 

All that said, the real message of FigurefTSl is how little difference any of this makes. The full 
range of source positions, once we are some distance downstream, is a few mm, and the total spread 
in the distance of the source to the MP, over all methods, is around 1%. 

Figure fT9l shows a practical application. The Burr Center neurosurgery beam is a single scattered 
beam using binary sets of Pb and Lexan slabs, manipulated by a computer, to scatter and range 
modulate the beam. As the slab configuration changes, so does the effective origin of protons, 
resulting in a small change in fluence. That change is compensated, to maintain a flat spread-out 
Bragg peak, by adjusting the in-beam time of each configuration. For Figure[T9l the equivalent 
source was computed according to the zbg (effective scattering point) method [8]. The others would 
give very similar source distances as we have just seen. 



4.11 Beam Fraction Contained in a Generalized Beam Ellipse 

The fraction of beam contained in the beam ellipse is of interest if, for instance, a degraded beam 
serves as the source for a magnetic beam line. It is obtained by integrating that part of the probability 
density Eq. (|20|) that lies within the ellipse. We can do that in closed form by transforming the ellipse 
into an erect ellipse (denoted with a prime) and then into a circle (double prime) , while the number 
of phase space points inside remains the same. 

Let us first define the ellipse more generally than we did in Eq.0Ol letting 

A x 2 - 2A Y x9 + A 2 9 2 = p 2 B (64) 



where p > is a dimensionless parameter. Following the methods of Sec. 14.31 this ellipse has 
bounding box x = py/~A~2 , 9 = py/Ao and area irp 2 

The effective extended source ellipse corresponding to the phase space probability density of 
Eq.[2D] is erect by definition. It has B' — B (emittance conserved in a drift) and A' Q = Aq (no 
scattering in a drift) but A' 2 = x 2 s = B 1 /A' = B/Aq, so in terms of the original As 

1 A x 2 + B6 2 /A 

P'(x,9)dxd9 = — X -=(T2 B dxd9 

27T V A B 



and the generalized beam ellipse is now erect with bounding box x' = p\J B j 'Aq , 6' = 9 = p\[A~Q and 
area Trp 2 y/~B. No phase space points cross the ellipse during a drift, so the beam fraction contained 
in the new ellipse is the same as before. 

Now define variables u, v such that the ellipse in u, v space is a circle: u = \/ Aqx, v = \J B / Aq 9. 
Let u 2 + v 2 = r 2 and replace du dv by r dr dip. Then 

P (r, </>) r dr d<p — = e 2Brdrd(f> (65) 

27TV B 
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The ellipse, now a circle of radius pV~B, still contains the same beam fraction and P" is still 
normalized cf. Eq. (|6]). Therefore the desired beam fraction is 

,2. rpVZ P"{pVB^) 
f P = J Q J o P M)rdr# = 1-e 2 = i _ _____ ( 66 ) 

The last equality says that if we take the ellipse as the contour where the probability density is 
down by e.g. 39.3% from its central value (corresponding to p = 1, the 'standard' ellipse), then 
the fraction of beam inside is also 39.3%. We have simply generalized to phase space a property of 
cylindrical Gaussians pointed out by Preston and Koehler [__] . 

Alternatively, one might wish to find the value of p and the corresponding emittance of an ellipse 
containing a stipulated fraction f p of the beam. From Eq. ([61)1 

P 2 = -21n(l-/,) (67) 

corresponding to emittance np 2 ^/B. 

For example, an ideal 230 MeV proton beam traverses a 10.7cm carbon degrader (1.8 g/cm 3 ), 
emerging with 150.2 MeV and emittance tt^/B = 19.63 7r mm mrad and entering a magnetic beam 
transport, whose designer defines the emittance ellipse as containing 95% of the beam. What is 
the emittance from her point of view? From Eq. (1571) . p 2 — 5.99 so the 2D emittance -Kp 2 \fB = 
118 mm mrad, six times the 'Fermi-Eyges' value. This happens in both transverse directions, and 
the energy distribution is also smeared out, so the transmission of a degrader/re-analyzer system is 
very small when the energy step is large. 



4.12 Summary 

The joint probability of a; and 9 (projected angle) at any depth z in a homogeneous slab 
has a Gaussian form fEq.l20"f involving three functions A n (z) which are moments of the 
scattering power T (Eqs.l2"Tl-l2l|. 

The combination B = A A 2 — A\ is positive if there is any scattering at all. 

The distribution of 9 irrespective of x is Gaussian with variance A$ (physical meaning 
of Aq, obvious from the definition of T). 

The distribution of x irrespective of 9 is Gaussian with variance Ai (physical meaning 
of A2, will be more obvious later). 

The rms spread of protons emerging from a slit at z about their mean angle is the angular 
confusion 9c — \J B/A2. It is independent of x (distance of the slit from the beam axis). 

The quantity conjugate to 9q is x c s = y^B/Ao , called x c g because it is the rms size of 
the effective extended source. 

The physical meaning of A\ is the covariance of x and 9 or < x 9 > . 

Ao, A2 and A\ (or B) can be regarded as the three parameters of the beam ellipse at 
z. The bounding box of the ellipse is given by v~ 2 and V~ o- The tilt of the ellipse is 
related to A\. Alternatively, the area enclosed by the ellipse (beam emittance) is ttVB 
(physical meaning of B). 

Therefore, given the three moments, the ellipse can be drawn (Eqs.l46l- [48|) . 

Applying the theory developed so far to a 127 MeV proton beam stopping in water, 
we find, at the maximum depth, an rms beam size \fA~2 that agrees with experiment 
(FigureEl). 

If a beam enters a finite homogeneous slab with moments A n its outgoing moments are 
A' n , given by Eqs. P91 - I5T]) . These equations can be iterated to transport an arbitrary 
(Gaussian) initial beam through mixed slabs. 

Eqs. (H9l - [5T|) imply that, in a drift, emittance is conserved (Liouville Theorem), as is Ao. 

Eqs. (|49l - [5~Tl) further imply that, in a thin scatterer, emittance increases by Eq. (|52|) 
(quadratic change = irx rms beam size x scatterer strength) and A2 is conserved. 
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Eqs. (|49l - l5T|) can be written in a differential form (Eqs.l54l-l56|) suitable for disordered 
heterogeneities (voxel-by- voxel treatment). 

It is sometimes useful to replace the entire beam line under consideration by an equivalent 
source drifting to the depth of interest in vacuo. Three such sources are defined: 

The effective extended source is that erect ellipse which drifts into the ellipse at the POL 
It exactly replaces the actual beam. 

The virtual point source is the point source from which the protons at the POI appear 
to emanate. 

The effective scattering point is the intercept of the back projected asymptote of o~ x (z) 
with the beam axis. 

The first two can be identified with intersections of the beam ellipse with its bounding 
box, while the last is the upper right hand corner of the box. Various equations can be 
written for parameters and positions of the sources. After any substantial drift in air or 
in vacuo the source distances from the measuring plane converge fFigurc lTS)) . 

The beam ellipse as defined conventionally (exponent in Eq.l2"0l equals —1/2) contains 
1 — exp(— 1/2) ~ 40% of the beam. If a different definition of emittance is required, 
Eqs. and (|57|) allow us to satisfy it. 

5 Beam Spreading in Matter 

As a simple but useful application of Fermi-Eyges theory we consider, in some detail, the transverse 
spreading of an ideal beam of protons or heavier ions in a homogeneous degrader. If the degrader 
happens to be water, that sets the limit on the smallness (as a function of depth) of targets that can 
be treated conformally in the ideal case. In practice, conformality will be worse because of non-ideal 
beams, heterogeneities, and non-elastic nuclear interactions. 

5.1 Theory of Preston and Koehler 

During the first decade of proton radiotherapy at the Harvard Cyclotron Laboratory (HCL), William 
Preston and Andreas M. Koehler studied, theoretically and experimentally, the limits imposed by 
multiple Coulomb scattering (MCS) on proton beams of very small initial cross section ('pencil 
beams' in current terminology). Unfortunately their manuscript |21j (hereinafter PK) was turned 
down[3 

In reading their paper one should realize that, though Fermi-Eyges theory already existed, they 
were unaware of it. Instead, they suggested an intuitive way of computing beam spreading, essentially 
deriving Eq. (|2"51) by direct physical reasoning. Using that model with the MCS theory of Bethe and 
Ashkin [23] they showed that the transverse spread at end-of-range of a proton beam stopping 
in any material is very nearly proportional to the range of the incident beam. The constant of 
proportionality depends on the material. They further showed that transverse beam spreading as a 
function of depth in any material obeys a universal law, for which they gave an analytic expression. 
That is, the basic shape of a x of the pencil beam vs. depth, expressed in appropriate reduced 
variables, is independent of the stopping material, the incident energy and (as we will show later) 
the properties of the stopping iono Let us re-derive these results in the language of Fermi-Eyges 
theory, closely following their reasoning. 

Figurel2"0l is a facsimile of PK's construction which makes Eq. (j2"3")l (recalling A 2 = < x 2 >) 
intuitively obvious. Consider the random kick, in a scatterer element dx, projected onto a measuring 
plane (MP) at S which can be anywhere downstream of the origin. A single kick has projected length 
(S — x)A9. The rms spread in the MP is the quadratic sum of such kicks. Passing from the sum 
to an integral over x, introducing the definition of scattering power, and changing to our present 
notation (x becomes z) we obtain Eq. (f23| . PK's results are all derived from this single equation. 

13 Years later Andy told me the journal (he could not remember which one) had deemed proton radiotherapy of 
insufficient general interest. The paper was a decade or two ahead of its time. 

14 A third important result (which we will not cover here) was their demonstration that, as the initial beam cross 
section is made smaller, the Bragg peak vanishes. The diminution of fluence on axis due to MCS eventually outweighs 
the increase in proton stopping power. In current language, transverse equilibrium is lost. 
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FIGURE I 



Figure 20: The construction of Preston and Koehler. The resultant displacement in the measuring 
plane equals the random sum of displacement vectors from infinitesimal elements of the slab. 



PK used a scattering power T derived from the MCS theory of Bethe and Ashkin [53] . We will 
instead use our 'differential Moliere' [TS] scattering power (Eqs. [TBl[T7|) P^l For a homogeneous slab, 
with z any depth < R\ , Eq. (|23j) then becomes 

a E 2 z 2 ~ [ z (z - z' f ^ , 

< x > z = — — fdu / —7 — 75— dz 



X S Jo (H 

(To avoid repetition we have anticipated heavy ions by writing a factor z\, the ion charge number, 
with every E s .) 

Though the scattering power is integrable numerically with fdM.{piv\, pv) included, to obtain 
results in closed form we have taken that slowly varying function outside the integral as an effective 
value fdM- in analogy to PK, let that be 

fdM = f<m(pivi, pv(z/2)) (69) 

that is to say, its value at half the depth of interest, call it /dM,z/2- Figure[2T1 shows the quality of 
that approximation. 

To integrate Eq. (|68[) we express the kinematic quantity pv in terms of z' using the excellent 
approximation^ 

(pv) 2 = a(pR) b = ap\R 1 -z') b = a p 1+k (R± — z') 1+k (70) 

where a and k depend on the stopping material, k varies inversely with scattering length Xs (or 
radiation length Xq) but is always small compared to 1 (Tabled]). Eq. (155]) becomes 

< x >z = ^7 /dM - /2 tf* I (r\ - zr>™ dz (71) 

If we desire the maximum (final) pencil beam size we set z = R\ and, integrating, obtain 



Elz\ 1 Rt k 



>Ri = — h /dM.Ri/2 TTT 7; T ( 72 ) 



15 Any T except Tpj^ yields almost identical results for water-like materials if only the transverse size of the beam 
at a given depth is considered |15| . 

16 Nowadays this is known as the 0veras approximation 1241 . Evidently PK were independently aware of it and the 
'weak 0veras' approximation used later. 
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depth (cm) 

Figure 21: Scattering power correction factor f^M for protons of five equally spaced ranges R\ (the 
largest corresponding to 250 MeV incident) stopping in water. The effective value /dM,z/2 used in 
evaluating Eq.[72] at each range is marked by a full circle. 
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4 He 2 



10 20 30 40 50 60 

pR (g/cm 2 ) 

Figure 22: Accuracy of pv reduction factor f\ = l./(z v /a) 1018 for three ions. Points are computed 
at 3, 6, 10, 20, 30 cm water and equivalent ranges (same proton energy) in Be, Al, Cu and Pb. 
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A slightly less good approximation is obtained by setting k = ('weak' 0veras approximation). In 
that case, pv should be refit as: 

(pv) 2 = ApiRt-z') (73) 



and Eq. ([72]) becomes 



E 2 z 2 

2 sir 
< X >Rj — — 7dM,Ri/2 



s £Ap 

This is PK's first result: a x at end-of-range for a given material is nearly proportional to range. If 
we evaluate A at Ri/2 and set \//dM,Ri/2 ~ 1-015 equal to 1 we are left with the very simple 



<Jx{Ri) 
Ri 



E s z\ 



2 0«)r 1 /2 V X s 



(75) 



which despite appearances is nearly independent of R\ . For protons in water with R\ = 5 cm it 
gives 2.37% , for 30 cm it gives 2.22%. The more complicated procedure based on Eq. (1721 (Tabled]) 
gives 2.25%. The three are equal for practical purposes. 

Continuing PK's derivation, in the weak 0veras approximation Eq. (|68[) becomes 



< x 2 > z 



X s 



/dM,z/2 ■ 



{z-z'f 



which can be integrated giving 



< x >, = 



Elz 2 
X S 



dM, z/2 



Ri_ 

Ap 



2 r 



(1 -i) 2 ln 



Ap J Ri 



+ t 



dz' 



1 - 1 



t-1 



Dividing Eq. ([77)1 by Eq. (|73j) and ignoring the difference between /, 



dM,z/2 



and /, 



dM,Ri/2 



(76) 

z/Ri (77) 
we find 



tr x (Ri) 



2(1- t) 2 In 



1 - t 



"I 1/2 



3r -2t 



z/Ri 



(78) 



This is PK's second result: pencil beam spreading is universal if a is normalized to its maximum 
value and depth is normalized to its maximum value (range). We shall find this remains true for 
any heavy ion. 

Table[T] lists parameters useful in evaluating the foregoing formulas for protons in 31 materials 
as well as the final constant of proportionality between a at end-of-range and range (see caption for 
details). From TableQ] 



protons in water: o^ru = 2.25% x Ri 
protons in aluminum: cr x , r x = 3.11% X Ri 



(PK Eq. (12c): 2.17%) 
(PK Eq. (12d): 3.18%) 



so our rederived numerical values are very near PK's (taking into account, of course, their use of 
a r = y2 x cr x rather than a x )- 



5.2 Generalization to Heavy Ions 

The results just obtained can be generalized to heavy ions. First we need to establish the range- 
energy relation for heavy ions and its consequence for pv. From stopping power theory [3] the mean 
projected range of any fast charged particle is 



R = 



i a r 0- 



c B Z 



L{fi) 



dr' 



(79) 



mc 2 is the particle's rest energy and z its charge number, t = E/mc 2 is the particle's reduced 
kinetic energy, eg is a universal constant while Z and A or (for compounds and mixtures) some 
combination of constituent Zs and ^4s, are atomic constants of the stopping material. L(j3) is the 



17 The lower limit of the integral is actually some very small reduced energy, a cutoff parameter whose exact value 
does not matter. 
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material 


P 


pX 


pX s 


k 


a 


rms 


.4 


rms 


rms 


<j x /Ri 




g/cm J 


g/cm 


g/cm 2 


% 




% 




% 


% 


% 


1 


2 


3 


4 




5 


6 




7 


8 


9 


10 


n 


HELIUM 


(0.166) 


94.32 


155. 


.55 


6.67 


3983 





.02 


4859 


3.85 


0.10 


1.67 


HYDROGEN 


(0.084) 


61.28 


145. 


.76 


6.11 


9221 





.01 


10552 


3.52 


0.02 


1.16 


BERYLLIUM 


1.848 


65.19 


92. 


.59 


6.98 


3285 





.04 


4093 


4.03 


0.14 


1.72 


POLYETHYLENE 


0.940 


44.64 


61. 


.79 


6.97 


4362 





.02 


5333 


4.02 


0.12 


1.89 


PILOTB 


1.020 


44.23 


59. 


.43 


7.01 


4097 





.03 


5036 


4.05 


0.14 


1.98 


POLYSTYRENE 


1.060 


43.72 


59. 


.15 


7.11 


3979 





.03 


4914 


4.11 


0.15 


2.01 


CARBLD 


1.800 


42.70 


56. 


.34 


7.18 


3582 





.05 


4463 


4.15 


0.18 


2.12 


CARBON 


1.900 


42.70 


56. 


.34 


7.18 


3582 





.05 


4463 


4.15 


0.18 


2.11 


LEXAN 


1.200 


41.46 


55. 


.05 


7.18 


3853 





.03 


4777 


4.15 


0.17 


2.10 


NYLON 


1.130 


41.84 


54. 


.73 


7.06 


4113 





.03 


5061 


4.08 


0.14 


2.05 


LUCITE 


1.190 


40.59 


53. 


.76 


7.18 


3943 





.03 


4881 


4.15 


0.16 


2.10 


KAPTON 


1.420 


40.56 


53. 


.21 


7.26 


3688 





.04 


4596 


4.20 


0.19 


2.16 


WATER 


1.000 


36.08 


46. 


.88 


7.29 


4007 





.04 


4970 


4.21 


0.18 


2.25 


AIR 


(1.205) 


36.66 


46. 


.76 


7.32 


3539 





.04 


4431 


4.23 


0.28 


3.05 


TEFLON 


2.200 


34.84 


43. 


.87 


7.49 


3305 





.05 


4178 


4.34 


0.26 


2.46 


QUARTZ 


2.640 


27.05 


32. 


.94 


7.89 


3247 





.06 


4156 


4.57 


0.35 


2.82 


MAGNESIUM 


1.740 


25.03 


30. 


.17 


7.98 


3161 





.06 


4065 


4.63 


0.40 


3.03 


ALUMINUM 


2.699 


24.01 


28. 


.75 


8.12 


3007 





.07 


3897 


4.70 


0.45 


3.11 


SILICON 


2.330 


21.82 


25. 


.93 


8.11 


3099 





.07 


4006 


4.70 


0.45 


3.25 


IRON 


7.874 


13.84 


15. 


.82 


8.85 


2576 





.13 


3453 


5.15 


0.86 


4.24 


COPPER 


8.960 


12.86 


14. 


.62 


9.13 


2446 





.13 


3321 


5.32 


1.01 


4.45 


NICKEL 


8.900 


12.68 


14. 


.42 


9.15 


2600 





.11 


3513 


5.32 


0.94 


4.36 


ZINC 


7.133 


12.43 


14. 


.09 


9.27 


2457 





.12 


3348 


5.40 


1.06 


4.56 


BRASS 


8.489 


12.30 


13. 


.89 


9.51 


2403 





.03 


3303 


5.46 


1.14 


4.58 


MOLYBDENUM 


10.220 


9.80 


10. 


.91 


9.77 


2168 





.15 


3033 


5.70 


1.54 


5.33 


TIN 


7.310 


8.82 


9. 


70 


10.13 


1995 





.16 


2843 


5.92 


1.96 


5.91 


GADOLINIUM 


7.900 


7.48 


8. 


04 


10.24 


1839 





.15 


2650 


5.97 


2.44 


6.69 


TANTALUM 


16.600 


6.82 


7. 


21 


11.18 


1717 





.20 


2566 


6.56 


3.04 


6.82 


MERCURY 


13.550 


6.44 


6. 


72 


11.23 


1628 





.21 


2450 


6.58 


3.44 


7.31 


LEAD 


11.350 


6.37 


G. 


63 


11.18 


1607 





.17 


2418 


6.54 


3.54 


7.49 


URANIUM 


18.950 


6.00 


6. 


12 


11.37 


1525 





.19 


2321 


6.65 


3.91 


7.71 



Table 1: Quantities related to the calculation of (J^^^/Ri in various materials. Eq.[75]is evaluated 
at five equally spaced ranges the largest corresponding to 250 MeV incident, and a straight line 
passing through the origin is fitted cf. PK Figure 3. Col. 11 is the slope of that line and col. 10 is 
the rms scatter of the five computed points about the line, a measure of the proportionality of cr Xj Rj 
to R\. Col. 2 is the density (parentheses: mg/cm 3 ), col. 3 is the mass radiation length and col. 4 is 
the mass scattering length. Cols. 5,6 are the Overas parameters from a two-point fit at i?i/3 and 
2Ri/3 and col. 7 is the rms scatter (at the five ranges) about that fit, a measure of the goodness of 
the approximation. Col. 8 is the 'weak Overas' parameter fitted at i?i/2 and col. 9, the rms scatter 
about that fit. 
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'stopping number', a function which, to a very good approximation, depends only on speed [3]. From 
Eq. (|10|) an ion's speed depends only on its r = E/mc 2 . Therefore, Eq. ([79| implies 



1 = 

as long as 



Rj (mc 2 /z 2 ) 



R p (mc? I ' z 2 ) 



t p = n or E p = Ei 

mi 



Introducing ai = mi/m p and noting that z p = 1 we find for any ion 



where R p is proton range as a function of energy.0 From this, if E p is proton energy as a function 
of range, 

2 

Ei{Ri) = aiE p (^R) (81) 

Given a proton range-energy function and its inverse (which must now extend to very large ranges 
because of z 2 /ai) we can now solve any range-energy problem for any ionF^I 

The 0veras approximations also generalize, as shown by Kanematsu [18) . The weak version 
Eq. (|73p reads, for protons, 

(Hp = A p (pR) p (82) 
Suppose a proton and an ion have the same speed v. From special relativity 



pc = mc 2 / \/l — p 2 
so momentum is proportional to mass and therefore 

(Hp = (Hi/ a i 

As already shown 

(pR) p = (z 2 /ai){pR)i 

Putting it all together we get the weak 0veras approximation for an ion in terms of the '0veras 
constant' of a proton 

(Hi 2 = (aiz 2 )(pR)iA p (83) 

Kanematsu |18) gives the analogous 'strong' approximation for water, also showing that it is much 
more accurate than the usual power law R = aE b . We will not need it here, however. 

To complete the generalization to heavy ions we need to consider the effect on scattering power. 
On elementary considerations (the Coulomb force) we must attach a factor zi to every occurrence 
of E s [H], which we have already done. Kanematsu's scattering power T^h [18] requires no further 
modification since its nonlocal correction factor depends only on the path integral of (1/Xq). 

We wish, however, to use our T^m (Eqs.[TH[l7]) whose correction factor /dM depends logarith- 
mically on the diminution of pv. f^M as it stands (Eq. ll7|) will not work since it was optimized for 
protons and characteristic (pv)s for ions of the same range are much larger, as we have just shown. 
Rather than opening Pandora's box and re-optimizing JdM for each ion let us apply an appropriate 
reduction factor fi to both p±v\ and pv letting 

/dM(piWi, pv) -> /divi(/i x (pivi), fi x (pv)) 

The form of /i is suggested by Eqs. (|82 [ I83 |) which, for protons and ions of the same range, give 

(Hp = (pv)i/(ziVai) (84) 



18 Of course R P (E) and E P (R) also depend on the stopping material. 

19 Our standard proton functions Range(energy, material) and Energy(rangc, material) assume kinetic energy in 
MeV and range in g/cm 2 , that is, mass range pR, despite consistently calling it 'range'. 
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A somewhat better approximation is 



fi 



1.018 



(85) 



Figurel22l shows it is good to approximately ±2% for materials, ions and ranges of clinical interest. 
That is good enough, because /dM only depends logarithmically on pv. For protons, fi = 1. 

Also involved in the derivation of T^m [E] was a condition on the characteristic large-angle 
(nuclear size) cutoff of Rutherford scattering. That lead to the simplification of X2/X1 an d thus to 
the introduction of scattering length X$- As easily shown, \2 < 1 is even better satisfied for heavy 
ions than for protons. 

Putting everything together, the derivations of the previous section go through exactly as before 
and we find that Eq. ([75]) for the R\ coefficient of rms beam spread and Eq. ([78]) for the universal 
beam spreading profile apply also to heavy ions. The influence of heavy ions in Eq. (|75p is found 
entirely in z\ and the fact that pv for heavy ions is much larger. 

Combining Eqs. (|84|) and (|75|) we also find, for ions and protons of the same range at any depth 



Kanematsu [22] obtains a.j , nearly the same. The reason for the disappearance or near 



disappearance of z is fundamental. Stopping power and scattering power both vary as z 2 . The extra 
Coulomb scattering is offset by the increased speed needed to get the same range. 

5.3 Experimental Tests 

A number of transverse scans of heavy ion beams vs. depth are available. Unfortunately they tend to 
be motivated entirely by clinical applications (beam commissioning), thus dominated by ion-optical 
effects or beam contamination [23] and not too useful as tests of the basic theory. After moderate 
effort we have found only the following three. 

5.3.1 Preston and Koehler 

The best (perhaps only) proton data remain those of PK [21 . Protons of approximately 158 MeV 
from the Harvard Cyclotron, degraded to lower energy as required and collimated to a 2 mm diameter 
beam, fell on a target of adjustable thickness. The beam profile exiting the target was measured 
using a very small remotely driven diode dosimeter |26] in conjunction with a plane parallel ion 
chamber beam monitor. The observed Gaussian distributions were fitted graphically to extract a. 
In-air a s at the same depths were subtracted in quadrature to correct for beam and detector size. 

PK's Figure 17 (our Figure l2"3"|) shows their universal curve along with data from two energies 
entering aluminum and one energy entering water. The data agree with theory (our Eqs. 1751 and [75|) 
within the experimental error of a few percent. 

5.3.2 Phillips 

This is a curious tale. A 1990 paper by Kraft [27] refers to 'data from M. Phillips, LBL' in a figure 
caption, with no further citation. Some of those data (for 1 Hi, 4 He2 and 12 Ce) were picked up by 
Hollmark et al. [28] who cited Kraft, and corrected for the size of the incident beam.0 Finally, 
Kanematsu [18] used the data referring to 'Phillips' measurements' and citing Hollmark. The data 
also show up in a Master of Science thesis by M. Granlund (2001) and perhaps elsewhere. 

LBL issued numerous technical notes at the time and it seemed unlikely that an experiment so 
difficult and fundamental would not have been written up or published. However, a library search 
turned up nothing. Aided by Prof. Kraft we located Prof. Phillips who, unfortunately, had discarded 
his documents from that period and remembered little of the episode: 

'. . . I suspect that the data that Gerhard Kraft originally published was a set of measure- 
ments I collated from the LBL particle therapy group (J. Lyman, W. Chu, T. Renner, B. 
Ludewigt) or actually just obtained from the treatment planning code data files (which 
of course would be data from the same group) . . .' 

20 That correction, obtained by extrapolation, corresponds to an unrealistically small beam diameter of 1.2 mm. It 
would, in any case, almost certainly have been done already by the original authors. 



&x,I = ^x., v l\fa\ 



(86) 
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Figure 23: Facsimile of PK Figure 17 whose caption reads "Dimensionless plot of the standard 
deviation due to scattering versus depth of penetration . . . Open triangles are experimental results 
for 112 MeV protons on aluminum; solid triangles for 158 MeV protons on aluminum; open circles 
for 127 MeV protons on water." 




He did, however, make us aware of the important LBL measurement to come next. 

The purpose of Kraft's paper was to promote heavy ion therapy and any reasonably correct data 
showing that heavy ions scatter less would have served. All things considered, we doubt that the 
data represent a measurement at all. However, in keeping with tradition, we compare the theory 
just derived with M. Phillips' data, read from Hollmark's graph, in Figure[24l We have not corrected 
for incident beam size. The agreement is quite good. 

5.3.3 Wong et al. 

This 29 is not a measurement of beam spreading but nevertheless extremely valuable as it confirms 
our understanding of MCS of heavy ions in an extreme case, 238 Ug2 scattering in Cu. 

Partly stripped U ions accelerated to 650 MeV/ai passed through a four element position-sensitive 
Si detector array, a 0.27 cm Cu target, and a second array. One assumes the U stripped completely 
on entering the Si. Each array made two x and two y measurements to establish incoming and 
outgoing track segments so the projected x and y deflections by the Cu were measured independently 
event- by- event. The final detector was thicker to provide a good AE measurement for background 
rejection. 

The analysis culminated in a Monte Carlo simulation using, in turn, three different MCS models 
namely 6*fr and ^Highland (in our terminology [TS]) and the full Moliere treatment of projected 
angle |13j . The MCS physics is thus somewhat buried, and the difficulty [TBI [15] of incorporating 
Highland's equation in a Monte Carlo was perhaps not fully appreciated. The authors concluded 
that Highland fit the final projected angle distribution best (it fits exactly), Moliere next (too broad) 
and #fr worst (even broader). Nevertheless 

'. . . all of the theories predicted a Gaussian distribution in satisfactory agreement with 
experiment, although the modified Fermi theory reproduced the data best . . .' 

as close as the authors come to giving an experimental error. 

Fortunately a much simpler analysis is possible using the information provided. Figure 6 gives 
target-in and target-out projected angular distributions for x and y. Scattering in Cu can be 
obtained by subtraction in quadrature, giving almost equal x and y angles which can be averaged 
and compared directly with various theoretical predictions using the mixed slab procedures and 
notation of [15]. In order of agreement with experiment (last column), with Kanematsu's 'differential 
Highland' scattering power [18] leading the pack: 



description 


mrad 


% 


^measured 


2.607 




#dH 


2.625 


0.66 


^Highland 


2.642 


1.35 


#dM 


2.683 


2.88 


^Hanson 


2.755 


5.67 



Of course, stopping powers or range-energy relations enter any such analysis or Monte Carlo. 
Those used by the authors were not directly described in the paper, and a citation proved unhelpful. 
Be that as it may, we used the heavy-ion formulas of the preceding section, the ICRU49 proton 
range-energy tables [3] , guessed at those Si detector thicknesses not given, combined each detector 
array and ignored air obtaining finally 



slab 


material 


cm 


MeV/a 


res cm Cu 


1 


Si 


0.3654 


650.0 


0.799 


2 


Cu 


0.2700 


587.2 


0.686 


3 


Si 


0.5324 


424.7 


0.416 


4 






310.5 


0.252 



Col. 3 is the slab thickness, Col. 4 is E entering the slab and Col. 4 is the residual range in Cu 
entering the slab. The only range-related benchmark in the paper is somewhat incidental, citing the 
maximum Cu that might have been used: 
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'. . . corresponding to uranium stopping in the final Y4 detector; it corresponds to a 
calculated copper target thickness of 3.5mm . . .' 

whereas we obtain 5.15 mm. We do not understand the discrepancy, but the good agreement obtained 
by Wong et al. suggest that their range-energy treatment was in fact correct. 

Our ranking of ^Highland and ^Hanson agrees with the authors' Monte-Carlo based conclusions. 
The amazingly good agreement of all four MCS theories with measurement must be fortuitous 
considering the likely experimental error. Still, the results of this experimental tour de force (2.6 mrad 
corresponds to 2.6mm in 1 m !) strongly reinforce our confidence in the MCS theory of heavy ions, 
as Wong et al. concluded. 

5.4 Summary 

We rederived PK's results [21] in Fermi-Eyges language. 

PK find the equivalent of Eq. (l23l) by reasoning that the rms transverse size of the beam 
at an arbitrary measuring plane is the quadratic sum of displacements, projected onto 
that plane, from infinitesimal elements of the stack. All their results follow from this 
single fragment of Fermi-Eyges theory. 

If the single scattering correction in the integral of Eq. ([23]) is replaced by its effective 
value and pv is replaced by the 0veras approximation then < a^ax > z — Ri) can 
be found in closed form fEq. [72|) . < x^ ax > is very nearly proportional to R\. 

If, instead, pv is replaced by the weak 0veras approximation, then < £^ ax > is exactly 
proportional to R\ and the end-point rms transverse beam size is proportional to range 
(PK's first result). The constant of proportionality is given by Eq. (f75|) . 

In the weak 0veras approximation the integral can also be evaluated at arbitrary depth 
leading (Eq. l78|) to a universal form for beam spreading in terms of normalized depth 
(PK's second result). It holds for protons of any incident energy traversing any material. 

The PK results can be generalized to arbitrary ions of mass/(proton mass) a\ and charge 
number Zj. Because stopping power depends only on these and speed, the range of an 
ion , given its kinetic energy, can be scaled from the proton range-energy relation (Eq. l80p 
and the energy of an ion, given its range, from the corresponding proton relation (Eq.l81|). 
The weak 0veras approximation scales similarly (Eq. 183]) . 

The scattering power for heavy ions depends on the single scattering correction. If, as 
in TdM, it is a function of pv it must be adjusted for the fact that pv for ions is much 
greater than for protons of the same range. Eq. [85] gives the appropriate factor to reduce 
pv to 'proton scale' and thus restore the validity of TdM ■ 

With all this done, the PK derivation goes through as before, extending their results to 
heavy ions. Beam spreading at any depth is less for a heavy ion than for a proton of the 
same range by l/^/ai (Eq. 186]) . 

There are few experimental tests of beam spreading in matter. PK validated their own 
theory with measurements |21] at several proton energies in aluminum and water. Data 
from Phillips at LBL have been widely cited but may or may not come from experiment. 
Wong et al. |29] is not a beam spreading measurement but nicely confirms multiple 
Coulomb scattering theory in an extreme case, uranium scattering in copper. 
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theory is not so well known. 
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A Analytical Geometry of the Ellipse 

Following Born and Wolf [30] we derive relations between parameters of the ellipse, centered on the 
origin, in several forms: the parametric form in the usual x,y frame, the parametric form in a u,v 
frame tilted to coincide with the principal axes, and the quadratic form in both frames. We also 
prove a formula for the area enclosed by a tilted ellipse. 




Figure 25: Ellipse geometry. 



A.l Tilted Ellipse 

The parametric form is 



By inspection this defines a closed curve depending on three parameters A, B, S , touching the boun- 
ding box 2 A x 2B at the points labeled in Figure l25l If 8 = 0, squaring and adding yields the 
quadratic form 

(5) 2 + (I)" - 1 m 



x = Acost (87) 
y = — Bs'm(t — S) 
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To express the more general tilted ellipse similarly we first rewrite Eqs. (|87 ll88l) as 

^ = cost (90) 

y 

— = — sin t cos 5 + cost sin 5 (91) 
B 

Taking [Eq. (|90)) x cos<5] and [Eq. (l90|) x sin 5 - Eq. (l9Tj) ] and simplifying we obtain 

— cos S = cos t cos <5 

x • x y -4. x 

— sin o — — = sin I cos o 
A B 

Squaring, adding and simplifying yields 

,/- 2 ™<i)(!) + (l) J = ™ 2s < 92 > 

the equation of a conic section, an ellipse or circle because the discriminant is negative: 
Asin 2 S/(A 2 B 2 )- A/{A 2 B 2 ) =- (2 cos 5 / AB) 2 < 

A. 2 Transformation to Principal Frame 

Let us now find the parameters of the same ellipse in the frame u, v tilted with respect to x, y by 
an angle ip such that u and v coincide with the principal axes of the ellipse (see Figure I25|) . The 
transformation from x, y to u, v is a counterclockwise rotation through ip 

u = x cos tp + y sin ip (93) 
v = — x sin ip + y cos ip (94) 

and the ellipse has the new and simpler form 

u = acos(t — to) (95) 
v = -bsm.(t-t ) (96) 

We assume a is the semimajor axis, a > b. to is not really an ellipse parameter since changing 
it does not change the shape or orientation of the ellipse. It is simply a shift in the origin of the 
time-like parameter t to ensure that Eqs. (f87ll88|) and Eqs. (j95ll96)) are capable of describing the 
same situation even with respect to the initial location on the ellipse. Since these equations are 
supposed to describe the same ellipse we have, on using Eqs. 



Expanding, we obtain 



acos(< — to) — Acost cos-0 — B sin(t — 6) smift 
-bsm(t — t ) = — Acost sin?/' — B sin (t — S) cos ip 



a (cost cos t + sin t sin t ) = ^4 cost cos tp — B sin tp (sin t cos S — cost sin S) 
b (sin t cos to — cost sin to) — ^4 cost sin -0 + Bcosip(sint cos 5 — cost sin <5) 

If these are to hold identically the coefficients of sin t and cos t must be equal: 

a cos to — A cos ip + B sin ip sin 5 (97) 

a sin to = — B sin ip cos S (98) 

b sin to = — A sin ip + B cos ip sin S (99) 

6 cos to = B cos ip cos 5 (100) 

These four equations contain all possible relations between a, 6, ip and A, B, S. One such is obtained 
by squaring each, adding and simplifying: 

a 2 + b 2 = A 2 + B 2 (101) 
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Another is obtained by taking Eq. ([97|xEq. ([T00])+Eq. ([98]lxEq. (|99| and simplifying: 

ab = ABcos6 (102) 
Finally, equating Eq. (|99|) /Eq. (|98|) to Eq. (|100|) /Eq. (|97|) . cross multiplying and simplifying we obtain 

2ABsin<5cos2V> = (A 2 - B 2 ) sin 2t/> (103) 
If we define an auxiliary angle a (where 0<a<7r/2)as 

t&na = B/A (104) 

we can rewrite Eq. (|103|) in the form 

tan 2 -0 = tan 2 a sin 5 (105) 

which gives the tilt angle ip in terms of A,B,S. Dividing Eq. (|102l ) by Eq. (|101[) and using the 
definition Eq. (|104l) we obtain 

2ab 2AB . n 

■ cos o = sin 2 a cos o 



a 2 + b 2 A 2 + B 2 
Defining a second auxiliary angle x (where < x < 7r/4) by 

tanx = V a (106) 

we can rewrite this as 

sin 2 x = sin 2 a cos <5 (107) 
A. 3 Summary and Explicit Procedures 

Eqs. ([571155)) and Eqs. (|55llM)l are the parametric form in the 'lab' and 'principal' frames respectively. 
Eqs. (|104[) and (| 1 06[) define useful auxiliary angles a and x, basically the aspect ratios in the two 
frames. Finally Eqs. (|101l) . (I105[) and (|107p relate a, b, ip to A, B, S using the auxiliary angles. Given 
A, B, S one may find a, b, ip using the procedure 

a = arctan£?/M. 

X = -5 arcsin (sin 2 a cos 6) 

•ip = .5arctan(tan2asin(5) 

A 2 + B 2 \ 1/2 



1 + tan 2 x 
b = a tan x 

Given a, 6, ^ one may find A, B, 6 using 

X = arctan6/a 



/tan 2 2V> + sin 2 2xV /2 

a = .5 arcsin — ■ 

V l + tan 2 2tp J 

S = arcsin (tan 2 ip/ tan 2 a) 

A _ c y /2 

V 1 + tan^ a J 
B = A tan a 

A. 4 Area Enclosed by the Ellipse 

The area enclosed 2 ^ by an erect ellipse of semiaxes a, b is easily obtained by integration: 



area. 



f 1 6(1 - {x/aff^dx = Aab [ (1 - u 2 f' 2 du 
o Jo 



,tt/2 

= Aab cos 2 Ode = nab (108) 



It is incorrect to say 'area of an ellipse'. The ellipse itself is a closed curve, having no area. 
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For a more general and less well known equation we note that, from Eq. (|92|) . B cos 8 is the y intercept 
of the ellipse and we already know that A is the maximum value of x. Therefore from Eqs. (|102j) 
and (TPS)) we find 

area — 7r x max — tt Xi n t Umax (109) 

the second version following from symmetry. Eq. (|109p reduces to Eq. (| 108[) for an erect ellipse. It 
is remarkable that Eq. (|109j) seems to be known only to accelerator physicists. We were unable to 
find it on the Web, unlike Eq. (j 108[) which is everywhere. 



List of Figures 

1 Top view of protons stopping in a water tank, with transverse scale exaggerated. 
Fermi- Eyges theory predicts the spatial and angular distributions of the full beam at 

any depth, and of the protons passing through the off-axis slit @ 

2 Model beam line for phase space analysis: three scatterers separated by voids or drifts, 

a collimator with a central slit, and a final drift 

3 Phase space evolution of the model beam (distributions at the entrance to each labeled 
element). For visual clarity the ellipse is drawn at 2.5 a rather than the conventional 
1 a; the standard ellipse would only contain about 390 of the 1000 phase space points. 

4 Phase space diagram of four protons acted on by a scatterer 

5 Phase space diagram of four protons acted on by a drift 

6 The beam ellipse, which encloses an area tt x[ nt a; ma x = tt ^max Xint 

7 Liouville's Theorem: the area enclosed by the beam ellipse (the beam emittance) is 
conserved in a drift 

8 Double scattered beam line with selected trajectories of protons that graze the down- 
stream edges of the patient collimator. SI represents one step of a range modulator. 

9 Phase space diagrams for the double scattered beam with x projections interleaved, 
(n.b. The setup used for these differs in two points from Figure[8l There is a collimator 
surrounding S2, and the 'patient collimator' is a half-beam block with one edge on 

axis rather than the symmetric one of Figured) I 

10 Projected angle distributions for 158.6 McV protons traversing 1cm of water. On 
a linear plot the Moliere/Fano distribution is indistinguishable from Gaussians using 
the Hanson or Highland 9 - On a log plot, the correct distribution peels away at 2.5ct, 

and is more than 100 x higher at 5<r LL2 

11 Stack of three homogeneous slabs for Fermi- Eyges discussion. Pb/Lexan/air simulates 
the upstream end of a double scattering beam line 

12 The beam ellipse in A n , B notation 

13 The beam ellipse in x, 9 notation 

14 Evolution of the beam ellipse of an ideal 127 MeV proton beam in near stopping 
thickness (0.97x range) water divided into five equal slabs. The dashed lines represent 
the uncertainty range of a measurement by Preston and Koehler |21j 

15 Evolution of the beam ellipse for a 230 MeV proton pencil beam entering a Pb/Lexan/air 
(respectively 0.234/14.6/85.1 cm) stack. The air is divided into 3 equal slabs 

16 Ellipses for the (top) effective extended and (bottom) virtual point sources 

17 Asymptotic 'effective scattering point' in a single slab 

18 Source locations v. depth of measuring plane using a 12 cm Lexan slab starting at 
z = 0, followed by air, and an ideal 160 MeV proton beam. Full circles use T^m ', 
dashed lines use Tic [S] ; open circles (zbg) use Highland's formula 8 (compare 
with z csp , full circles of same size) [25 

19 Practical application of the effective scattering point. The first scatterer of the Burr 
Center 'STAR' beam uses Pb and Lexan degraders controlled by a computer to pro- 
duce the desired SOBP. As the degraders change so does the effective scattering point 
(filled squares). Note the jump at the major binary transition. The related change in 
1/r 2 is compensated by adjusting the dwell time of each configuration in the beam. . [25 

20 The construction of Preston and Koehler. The resultant displacement in the measur- 
ing plane equals the random sum of displacement vectors from infinitesimal elements 

of the slab EH 
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21 Scattering power correction factor fdM for protons of five equally spaced ranges R\ 
(the largest corresponding to 250 MeV incident) stopping in water. The effective value 
/dM,z/2 use d in evaluating Eq.[75]at each range is marked by a full circle [3C 

22 Accuracy of pv reduction factor j\ = 1./ (zy/a) 1 - 018 for three ions. Points are computed 
at 3, 6, 10, 20, 30 cm water and equivalent ranges (same proton energy) in Be, Al, Cu 

and Pb E 

23 Facsimile of PK Figure 17 whose caption reads "Dimensionless plot of the standard 
deviation due to scattering versus depth of penetration . . . Open triangles are experi- 
mental results for 112 MeV protons on aluminum; solid triangles for 158 MeV protons 



on aluminum; open circles for 127 MeV protons on water." |35 

24 Comparison between 'data from M. Phillips, LBL' (see text) and PK theory general- 
ized to heavy ions (Eqs.[7Sl[75]) (lines) 

25 Ellipse geometry 
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